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ABSTRACT 

Observations of SN 1987A by the Chandra High Energy Transmission Grat- 
ing (HETG) in 1999 and the XMM-Newton Reflection Grating Spectrometer 
(RGS) in 2003 show very broad (v-b) lines with a full-width at half-maximum 
(FWHM) of order 10^ kms~^; at these times the blast wave was primarily in- 
teracting with the Hll region around the progenitor. Since then, the X-ray 
emission has been increasingly dominated by narrower components as the blast 
wave encounters dense equatorial ring (ER) material. Even so, continuing v-b 
emission is seen in the grating spectra suggesting that interaction with Hll re- 
gion material is on-going. Based on the deep HETG 2007 and 2011 data sets, 
and confirmed by RGS and other HETG observations, the v-b component has a 
width of 9300 ±2000 kms"^ FWHM and contributes of order 20 % of the current 
0.5-2 keV flux. Guided by this result, SN1987A's X-ray spectra are modeled as 
the weighted sum of the non-equilibrium-ionization (NEI) emission from two sim- 
ple ID hydrodynamic simulations, this "2xlD" model reproduces the observed 
radii, light curves, and spectra with a minimum of free parameters. The inter- 
action with the Hll region (pinit ~ 130 amucm"^, ±15 degrees opening angle) 
produces the very-broad emission lines and most of the 3-10 keV flux. Our ER 
hydrodynamics, admittedly a crude approximation to the multi-D reality, gives 
ER densities of ~ 10^ amucm"^, requires dense clumps (x5.5 density enhance- 
ment in ~ 30 % of the volume), and it predicts that the 0.5-2 keV flux will drop 
at a rate of ~ 17 % per year once no new dense ER material is being shocked. 
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troscopic — X-rays: general 
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INTRODUCTION 



Since its explosion as the first optical SN observed in the year 1987, SN 1987A has 
been the subject of intense study at every possible wavelength. Its locati on in the Larg e 
Magellanic Cloud (LMC), at a relatively close distance of 51.4 ±1.2 kpc (jPanagial Il999l ). 
has enabled us to study details regarding the aftermath of a SN explosion not previously 
possible; we have learned and will continue to learn more about the physics of SNe and the 
formation of supernova remnants (SNRs). 

The X-ray emission from SN 1987A has likewise been studied with every available satel- 
lite. However, SN 1987A is by far the intrinsically dimmest X-ray SN to be observed, at 
least in its f irst decade, and is only now a pproaching the luminosity of other comparable-age 
X-ray SNe (iDwarkadas &: Gruszkdl2012l ). It is clear that the only reason it was detected in 



early observations at high X-ray energies (6-28 keV) with the Ginga satellite fllnoue et al. 



19911 ) was because of its proximity. The emission appeared to decrease over the next couple 
of years, and it appeared as though SN 1987A would behave like every other X-ray SN then 
known, with an X-ray emission decreasing with time. 



Just over 3 years after its expl osion however, the X- ray (IHasinger et al.lll996l . 0.5-2 keV 
ROSAT observations) and radio ( iGaensler et al.l 119971 . 1.4-8.6 GHz with the Australian 
Telescope Compact Array) emission from SN 1987A bega n instead t o inc rease with time. 
The increase in radio emission was initially interpreted by IChevalierl ( 119921 ) as arising from 
the interaction of the SN forward shock (FS) with the wind termination shock of the blue 
supergiant (BSC) wind of the progenitor. However this would not result in a continued 
increase but would instead be limited in time. The observations indicated otherwise, and 
the X-ray emission continued to increase. The subsequent increases have suggested that the 
interaction is with a much higher density and mo re extended region than expecte d from the 
wind termination shock. This was interpreted by IChevalier fc Dwarkadad ( 119951 ) as ionized 
red-supergiant (RSG) wind emitted during a pre-BSG phase, and subsequently swept-up by 
the fast BSC wind. The ejecta interaction with the resulting Hll region produces the usual 
two-shock structure with the FS moving into the Hll region and a reverse shock (RS) moving 
back into the ejecta; a contact- discontinuity (CD) d efines the ejecta-HlI bou ndary. Such a 
model was used to explain the early X-ray emission (IBorkowski et al.lll997al ) and predicted 
line widths of 5000 kms~^ or more. 

The X-ray emission since then has continued to rise, as shown in Figure 1 for the 0.5- 
2 keV and 3 -10 k eV bands. A review of the X-ray emission over the first 20 years is given 



m 



McCrayl (120071 ). The hard X-ray and radio emission have continued to rise together. 



but after about 6000 days the soft X-ray light curve has steepened still further, indicating 
a source of emission that leads predominantly to soft X-rays. This has been interpreted 
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by iMcCrayl (120071 ) as due to the FS interacting with dense fingers of emission arising from 
the inner edge of the equatorial ring (ER) surrounding SN 1987A. From days 6500 to 8000 
the soft X-ray fiux continued to grow at a roughly linear rat^, that is dF/dt constant. 
Most recently, iPark et al.l (120111 ) have reported a fiattening of the X-ray light curve since 
day ~ 8000, indicated in Figure [1], and have suggested the possibility that the FS has now 
propagated beyond the majority of the dense inner ring material. This could lead to further 
changes in the light curve and, hence, future measurements are eagerly anticipated. 

Even considering its low luminosity, the proximity of SN 1987A has allowed for high reso- 
lution grating observations to be taken. These observations, including recent ones from 2011, 
are summarized and analyzed in §[2l There is much to learn from this set of high-resolution 
X-ray spectra; however, here we focus on the existence and evolution of a "very-broad" (v- 
b) emission line component in the data, § 12.41 This v-b component is pe rsistent over the 



l ast dec ade and it is likely the continuation of the br oad lines expe c ted by iBorkowski et al. 



(Il997al ) and measured in the Chandra HETG data of iMichael et al.l (120021 ). In §[3] we review 
and present our hydrodynamic modeling of SN 1987A which consists of two contributions: 
the Hll material above and below the equatorial plane and the dense equatorial ring ma- 
terial. We compare the results of the model with observed quantities in §IH showing that 
the v-b component naturally arises from the shocked Hll material and is the main source of 
the 3-10 keV fiux. Also, our very simplified hydrodynamics of the ring interaction is used 
to show how the 0.5-2 keV light curve would behave if the FS has indeed gone beyond the 
densest ring material. In the final section we summarize our conclusions. 



^This was first pointed out by D.N. Burrows in discussions during the preparation of iPark et al.l ( 2011 ) 
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Fig. 1. — SN 1987A X-ray light curves with grating observation epochs indicated. The 0.5- 
2 keV data (black, red) show an initial exponential behavior whose slope increases after 
day ~4500 and then decreases from day ~6300, possibly leveling off at the recent epoch 
~8800 days. In constrast, the 3-10 keV flux (green) has more closely followed a consistent 
exponential behavior to the present. For reference, the sum of an exponential (beginning at 
day 2000) and a linear ramp increase (from day 6000 to day 8200), is included (magenta, 
blue) to hi ghlight the ge r ieral p hases of the 0.5-2 keV light c urve. The ROSAT fluxes are 
taken from iHaberl et al.l (120061 ) and the ACIS values are from lPark et al.l ( 120111 ). 
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GRATING OBSERVATIONS AND DATA REDUCTION 



The grating observations of SN 1987A are listed in Table [H including those from in- 
struments on Chandra (High- and Low- Energy Transmission Gratings: HETG, LETG) 
and XMM-Newton (Reflection Grating Spectrometer, RGS). The observation epochs are 
marked o n the SN 1987A ligh t curve of Figure [H showing that only the HETG-99 ob- 
servation (IMichael et al.l |2002| ) was taken before day 5500 after which the flux from the 
shock-protrusion collisions became significant . Along with HE TG-99, the deep non-grating 
Chandra ACIS-00 observation at 5036 days (iPark et al.l |2002| . Obsid 1967) gives the best 
pre-coUision X-ray spectra (albeit at only CCD resoluti on) and is included in the table as 
well. With the coll i sion g etting underway, the RGS-03 jHaberl et al.ll2006h and LETG-04 
( iZhekov et al.ll2005l . |2006| ) observations cover a transition period. Since 2007, with most of 
the X-ray flux due to the shock -ring interaction, there have been regu lar grating observa- 
tions made with b oth Chandra ( jPewey et al.l l2008l : IZhekov et al.l 120091 ) and XMM-Newton 
dsturm et allboioh . 



For our data analysis, we use standard spectral extractions from the observations as 
described i n ^12.11 fc ^12. 2[ The vario us analyses, described in § 12.31 & § 12. 4[ are then carried 
out in ISIS (IHouck fc Denicolall2000l ) using a common set of scripts that include some simple 
tailoring of the energy ranges and spatial-spectral parameters (see Appendix \^ based on 
the type of data, e.g., HETG vs. RGS. 



2.1. Chandra Data 



Public Chandra grating data and their spectral extractio ns are available and were ob- 
tained from the TGCat archive^ (jHuenemoerder et al. 2011 ). All of these had been re- 
processed in 2010 or later and therefore the downloaded pha, arf and rmf files could be 
used without any further re-processing. Multiple Obsld's at the same epoch were combined 
in ISIS by summing the pha counts, creating an exposure- weighted arf, and summing the 
exposure times. To reduce background from the less sensitive front-illuminated CCDs, the 
pha and arf for wavelengths longward of 17.974 A were set to zero in the plus [minus] order 
of the MEG [LEG] data. The plus and minus orders are then combined, even though they 
do have sm all but significant differences in their li ne- widths because of the resolved nature of 
SN1987A (IZhekov et al.ll2005l : iDewey et al.ll2008l ). The composite line-shapes are accounted 
for in the model by a smoothing function that has parameters based on the observed line 



2 



http : //tgcat . mit . edu/ 
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widths, § 12.31 and Appendix |X1 In this way we reasonably approximate the combined-orders' 
hne shapes. Thus, further analysis is carried out on a single spectrum in the case of LETG 
data and on two spectra (MEG and HEG) for HETG observations. 

Some HETG observations are indicated with a lower-case "hetg-YY" designation in the 
table. These spectra are from the combined "monitoring observations" (PI - David Burrows) 
in the year 20YY. Starting in 20 08 the HETG was inserted for SN 1987A observations in order 
to reduce the effects of pileup (jPark et al.ll201ll ): as a by-product, additional HETG high- 
resolution data are taken. However, to obtain the shortest exposure times only a subframe 
of the A CIS is read out; because of this the outer portions of the HETG dispersed "X" 



pattern ( ICanizares et al.l l2000l |2005[ ) fall outside of the subframe, removing response at 



longer wavelengths. Instead of centering SN 1987A in the subframe, it is offset in cross- 
dispersion to get more of the sensitive MEG minus order in the subframe and this extends 
the response to just include the OVIII line at ~0.65 keV. 

The most recent HETG observations of SN 1987A, HETG-11 in Tabled], were carried 
out during 2011 March 1-13 as part of the GTO program (PI - Claude Canizares). The 
proprietary data were retrieved from the Chandra archive and scripts from TGCav^ were 
used to carry out the spectral extraction using standard CIAO tools; specifically CALDB 
4.4.5 and CIAO 4.3 were used to process the data. 

The non-grating ACIS-00 observation, Obsid 1967, was downloaded from the archive 
and processed through level 1 event filtering with CIAO tools. The psextract routine was 
then used to generate pha, arf, and rmf files for a ~6 arc second diameter region around 
SN 1987A; a background extraction was made but not used because it had very few counts. 



2.2. XMM-Newton RGS Data 



The European observatory, XMM-Newton, has been use d to study SN 1987A a t high 
spectral resolution via its reflection-grating spectrometer, RGS (Iden Herder et al.ll200ll ). The 
data and processir ig steps for the four epochs of RGS data taken through January 2009 have 
been presented in ISturm et al.l (120101 ). Two subsequent observations, RGS-10 and RGS-11, 
have been carried out and were similarly processed; XMM-Newton SAS 10.0.0 was used 
in (re-)processing all RGS data. When read into ISIS the RGS-1 and RGS-2 counts (and 
backgrounds) are kept separate and jointly fit in further analysis. 



^ See Help=^> Software on the TGCat web page: 



http : //tgcat . mit . edu/ 
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Table 1. High-Resolution X-ray Observations of SN1987A to an Age of 25 Years 





Epoch 


Obsld(s) 


Expos. 


SN Age 


F(0.47-0.79) 


F{0.79-1.1) 


F(l. 1-2.1) 


Data ID " 


(year) 




(ks) 


(days)'^ 


(Cf 


;s)'= 


(cgs) " 


(cgs) " 


HETG-99 


1999 76 


124 1387 


116.2 


4606 


0.94 


[12%] 


0.57 


[9%] 


0.85 


[5%] 


ACIS-00'* 


2000.93 


1967 


98.8 


5036 


0.67 


[2%] 


0.71 


[2%] 


1.01 


[2%] 


RGS-03 


2003.36 


0144530101 


111.9 


5920 


4.2 




3.6 




3.4 




LETG-04 


2004.66 


4640, 5362, 4641, 6099, 5363 


288.8 


6396 


5.4 


[2%] 


6.2 


[2%] 


5.2 


[1%] 


RGS-07 


2007.05 


0406840301 


110.5 


7268 


10.0 




13.8 




12.2 




HETG-07 


2007.23 


8523, 8537, 7588, 8538, 7589, 
8539, 8542, 8487, 8543, 8544, 
8488, 8545, 8546, 7590 


354.9 


7335 


13.4 


[3%] 


15.5 


[1%] 


12.6 


[1%] 


LETG-07 


2007.69 


9581, 9582, 9580, 7620, 7621, 
9591, 9592, 9589, 9590 


285.2 


7503 


13.4 


[2%] 


18.2 


[1%] 


15.4 


[1%] 


RGS-08 


2008.03 


0506220101 


113.7 


7627 


12.7 




18.0 




14.5 




hetg-08 


2008.50 


9144 


42.0 


7799 






19.6 


[3%] 


18.3 


[2%] 


RGS-09 


2009.08 


0556350101 


101.4 


8012 


15.3 




22.2 




17.9 




hetg-09 


2009.30 


10221, 10852, 10853, 10854, 
10855; 10222, 10926 


129.7 


8091 






23.2 


[2%] 


21.2 


[1%] 


RGS-10 


2009.93 


0601200101 


91.6 


8320 


16.4 




25.2 




21.7 




hctg-10 


2010.46 


12125, 12126, 11090; 
13131, 11091 


118.3 


8515 






25.0 


[2%] 


24.8 


[1%] 


RGB- 11 


2010.93 


0650420101 


65.7 


8687 


17.4 




29.0 




24.9 




HETG-11 


2011.17 


12145, 13238, 13239, 12146 


177.7 


8774 


15.8 


[3%] 


27.6 


[1%] 


27.4 


[1%] 


hctg-11 


2011.47 


12539; 12540, 14344 


101.3 


8883 






25.2 


[2%] 


27.1 


[1%] 


RGS-12 


--'11.92 


0671080101 


-68 


~9048 














hetg-12 


-^'12.22 


13735, 14417 


-75 


~9157 















''The lower-case "hetg" data sets indicate ACIS monitoring observations where a narrow readout is used; this reduces the 
wavelength coverage of the HETG. 

''Epoch gives the observation-averaged date of the data set; SN Age is in days-since-explosion given by 365.25*( Epoch — 
1987.148). 

F{Ei-E2) is the observed fiux from Ei keV to E2 keV in units of 10"^'' erg/s/cm^. It is model-independent for the Chandra 
data sets with the statistical error given in brackets. The RGS fluxes are based on the 3-shock model fits (ii l2.3l l: statistical 
errors are not given since model-fitting errors of 5-10% dominate. 

''This non-grating ACIS observation is included here for reference because it is the deepest very-early observation. 
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2.3. Fluxes and Three-Shock Model Fitting 

The Chandra grating observations can robustly determine the observed flux (and sta- 
tistical error) in energy bands directly from the data without using a model. Instead of the 
the usual coarse "flux in the 0.5-2 keV band," fluxes in Table [1] are given for three bands 
covering this range. The band boundaries are chosen to avoid strong lines and are dominated 
by lines of N&O, Fe&Ne, and Mg&Si, respectively. Because of their reduced wavelength 
coverage the "hetg-YY" data sets do not allow an accurate flux measurement in the lowest- 
energy band. For the RGS data sets the fluxes given in the table are from 3-shock model 
flts (described below) because the RGS spectral gaps and background counts complicate a 
direct-from-the-data approach. 

Comparing HETG-11 fluxes to the previous HETG-07 ones shows that the main devi- 
ation from the trend of continued flux increase has happened primarily in the lowest energy 
range: while the 0.79-1.1 and 1.1-2.1 keV fluxes have increased by factors of 1.78 and 2.17, 
respectively, the 0.47-0.79 keV flux has increased by a smaller factor of only 1.180 Likewise, 
similar spectral changes are seen in the recent RGS data of 2010 December ("RGS-11"), 
continuing trends seen earlier; for example, whereas the Ne X flux increased by 36 % from 
RGS-08 to RGS-09 the lower-energy NVII flux increased by only 10% in the same period 



(ISturm et al.ll2010f l. 



A nominal non-equilibrium ionization (NEl) model was flt independently to each epoch's 
data with the goal of capturing all of the relevant lines in that epoch's spectrum. In order 



to allow for the wide range of shock temperatures seen in SN 1987A (iZhekov et al.l 12009 
Figure 2), the model consists of the sum of three plane-parallel shock components with 
common ( "tied" ) abundances. Physically, a single plane-parallel shock model represents the 
emission from a uniform slab of mat erial into which a shock wave has propagated for some 



finite time, tg ( iBorkowski et al.ll200ll ). The shocked portion of the slab has uniform electron 
temperature, kT^, and density, rig, but has a linear variation in the time since shock passage, 
going from ~ immediately postshock to a maximum of tg for the earliest-shocked material. 
Thus, the plane-parallel shock emission can be parameterized by a normalization (oc nennV), 
the abundances, the postshock electron temperature, and the maximum ionization agelfl. 



The apparent low-energy decrease may have a contribution due to cahbration issues, in particular with 
the ACIS contamination modeling. An HETG observation of the SNR E0102 was taken on 2011 February 
11, just a month before the SN 1987A data; its initial analysis shows calibration changes in the oxygen lines, 
but only of order 20%. 

^ Specifically, the maximum ionization age is specified by the Tau_u parameter of the vpshock model in 
XSPEC; there is also a Tau_l parameter which in keeping with our physical description is set to 0. 
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The sum of the three shock components is then smoothed with a Gaussian function 
to account for hne-broadening and the finite size of SN 1987A. Finally, there is an overall 
photo-electric absorption term to account for Galactic and LMC absorption. The model is 
implemented in ISIS, making use of the XSPEC spectral models library; the model specifica- 
tion (Expression lAip and a description of the smoothing par ameters are given in Appendix A. 
Many of the model parameters are fixed as in previous work (jZhekov et al.ll2009l ): the column 



densityl^ {Njj = 1.3x10^^ cm ^), the redshifts (286.5 kms ^), and the abundances of elements 
without strong, visible lines (H=l, He=2.57, C=0.09 N=0.56, Ar=0.54 Ca=0 34, Ni=0.62); 
as previously the solar photospheric abundances of lAnders fc Grevessd (119891 Table 2) are 
used as the reference set. There remain the free parameters consisting of 3 normalizations, 
3 temperatures, 3 ionization ages, and the abundances of O, Ne, Mg, Si, S, and Fe. 

The parameter set was reduced further by fixing the middle temperature at a value near 
the (logarithrnic) ce nter of the emission measure distribution of SN 1987A, kT^i^ = 1.15 keV 
(jZhekov et al.ll2009l . Figure 2). Since the ~1 keV emission had grown from 2004 to 2007, 
the expectation, borne out in Tables is that this component would have a growing 

contribution to the model, especially at later epochs. To test the (in)sensitivity of the fit 
values to the exact value of fcTmid, the HETG-11 data were fit with additional kT^i^ values 
of 1.0 & 1.3 keV. All of the HETG-11 fit parameters showed small variations over the test 
set of mid-temperature values, {1.0, 1.15, 1.3} keV. Some examples are: kTio takes values 
{0.50, 0.54, 0.56}, the fcThi norm varies as {2.4, 2.0, 1.5}, the silicon abundance varies as 
{0.41, 0.39, 0.38}, and the changed slightly as {1069, 1088, 1111} (for 560 degrees of 
freedom.) These small, gradual variations demonstrate that the exact choice for the middle 
temperature is somewhat arbitrary and nearly degenerate with other parameters. 

Two more p arameters are r emove d by constraining the ionization ages to cover a factor of 
2, e.g., as seen in lZhekov et al.l (|2009[ ). specifically setting: rio = -\/2 Tmid and = Tmid / a/2. 
Hence, in addition to the abundances there are 6 free parameters: kT\o, r^i^, kT^^i, and the 
3 normalizations, A'lo, A'^i.is, A'^hi- This is the same number of parameters as for a general 
2-shock model (2 norms, 2 fcT's, 2 r's) and gives slightly better fits. 



^ The column density could have a measureable local component that would change as the SN 1987A 
system evolves: taking a density of 1x10"* Hcm^'^ with a 0.01 pc path length sets an A^h local scale of 
'^0.3x10^-'^. Clumping, the ionization state, and the geometry of SN1987A have to be taken into account 
as well. In terms of data analysis, the Nh is very degenerate with the low-energy line abundances (N, 
O) and with the continu um teperature s and norms. For our purposes we fixed Nh at the two-shock- 
model value determined in IZhekovet all|2009l) . roughly half of which is due to the Galactic contribution of 
-0.6xl02icm-2. 
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Fig. 2. — HETG-11 data and model. The data (black) are globally fit by a 3-vpshock 
model (red). Removing the high-temperature component gives the blue curve and just the 
low-temperature component is shown in green. 



The 3-shock model was fit to the HETG and LETG data in the 0.6 to 5 keV range, 
giving good overall fits, Figure [2l For fitting, the spectra were binned to a minumum of 
20-30 counts per bin and a minimum bin width of ~ 0.05 A. The HEG data are ignored 
below 0.78 keV (>15.9A) and ignored completely for the HETG-99 data. Because of their 
reduced wavelength coverage, several changes were made for the "hetg-YY" data sets: they 
were "noticed" in the range 0.69-5 keV for the MEG (1.25-5 keV for the HEG), the O 
abundance was fixed, and the kT\o value was fixed (these latter two were set to values based 
on the HETG-07 and HETG-11 fits.) The RGS data sets do a good job of covering the N 
line ~0.5 keV but have reduced sensitivity above the Si XIV line, so these data were fit in 
the range 0.47 to 2.20 keV with the N abundance free, the S abundance frozen at a nominal 
value (0.36), and the kT^ value fixed. The ACIS-00 spectrum was fit over the 0.4 to 8.1 keV 
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range with N free. 

The values for the 3-shock fit parameters and their 1-a confidence ranges are given in 
Tabled As expected from their low fluxes, the two pre-collision data sets (HETG-99, ACIS- 
00) show comparatively low normalizations. Note that the RGS values and the HETG/LETG 
values differ significantly at similar epochs; this is probably due to near- degeneracies in the 
model which amplify the differences between the spectrometers in terms of their energy-range 
coverage, statistical weighting (counts vs. energy), and line-spread functions. However they 
do show the same general trends, e.g., toward larger Tmid values at more recent epochs. In 
the next section we use these 3-shock fits as good starting points to look for and measure 
the presence of the v-b component. 



Table 2. Parameters of the Global 3-Shock Model Fits 







Normalizations 


a 




T mid 






— Relative abundances ^ — 




Data ID 




A^i.is 




(keV) 


(10"^) 


(keV) 


N 


O 


Ne 


Si 


Fe 


HETG-99 


0.17 ±35% 


[0] 


0.13 ±12% 


[0.74] 


2.1 ±342% 


[4.3] 


[0.56] 


0.12 ±286% 


0.13 ±185% 


0.13 ±135% 


0.10 ±160% 


ACIS-00 


0.18 ±20% 


[0] 


0.12 ±9% 


0.73 ±11% 


8.2 ±19% 


4.2 ±14% 


1.6 ±63% 


0.47 ±50% 


0.61 ±42% 


0.68 ±27% 


0.12 ±45% 


RGS-03 


3.6 ±12% 


[0] 


0.14 ±5289% 


[0.50] 


4.1 ±46% 


[4.3] 


0.24 ±41% 


0.053 ±34% 


0.11 ±21% 


0.42 ±184% 


0.027 ±23% 


RGS-07 


3.0 ±6% 


0.35 ±36% 


1.1 ±11% 


0.51 ±2% 


4.8 ±11% 


[4.3] 


1.1 ±12% 


0.15 ±5% 


0.47 ±9% 


0.61 ±19% 


0.28 ±7% 


RGS-08 


3.4 ±3% 


0.72 ±15% 


1.3 ±2% 


0.52 ±2% 


4.7 ±1% 


[2.7] 


1.1 ±4% 


0.15 ±4% 


0.40 ±7% 


0.37 ±21% 


0.30 ±5% 


RGS-09 


5.3 ±3% 


1.4 ±9% 


1.2 ±12% 


0.53 ±1% 


7.9 ±6% 


[2.7] 


1.2 ±0% 


0.17 ±0% 


0.44 ±5% 


0.35 ±19% 


0.25 ±3% 


RGS-10 


5.1 ±3% 


1.6 ±7% 


1.5 ±0% 


0.54 ±1% 


10 ±0% 


[2.7] 


1.5 ±0% 


0.20 ±0% 


0.53 ±5% 


0.57 ±7% 


0.29 ±2% 


RGS-11 


4.7 ±3% 


1.3 ±0% 


2.1 ±6% 


0.57 ±1% 


15 ±0% 


[2.7] 


2.0 ±0% 


0.31 ±0% 


0.76 ±0% 


0.61 ±16% 


0.39 ±2% 


LETG-04 


2.3 ±4% 


0.64 ±10% 


0.28 ±8% 


0.48 ±2% 


3.7 ±8% 


5.3 ±9% 


[0.56] 


0.090 ±9% 


0.26 ±6% 


0.30 ±6% 


0.13 ±6% 


HETG-07 


5.2 ±1% 


1.7 ±2% 


0.64 ±2% 


0.47 ±0% 


3.4 ±0% 


4.3 ±1% 


[0.56] 


0.096 ±0% 


0.22 ±0% 


0.31 ±0% 


0.16 ±0% 


LETG-07 


5.2 ±2% 


1.4 ±2% 


1.2 ±2% 


0.48 ±1% 


3.7 ±4% 


2.4 ±1% 


[0.56] 


0.090 ±3% 


0.27 ±2% 


0.34 ±1% 


0.19 ±2% 


hetg-08 


6.6 ±7% 


1.6 ±37% 


1.5 ±31% 


[0.52] 


3.6 ±18% 


2.4 ±18% 


[0.56] 


[0.11] 


0.20 ±10% 


0.32 ±8% 


0.14 ±10% 


hetg-09 


5.5 ±6% 


2.2 ±10% 


1.6 ±9% 


[0.52] 


5.9 ±12% 


2.7 ±6% 


[0.56] 


[0.11] 


0.30 ±7% 


0.37 ±5% 


0.22 ±5% 


hetg-10 


5.9 ±7% 


2.8 ±14% 


1.9 ±16% 


[0.52] 


5.7 ±13% 


2.7 ±12% 


[0.56] 


[0.11] 


0.27 ±9% 


0.37 ±5% 


0.19 ±5% 


HETG-11 


6.2 ±1% 


3.1 ±0% 


2.0 ±0% 


0.54 ±1% 


7.0 ±0% 


2.7 ±0% 


[0.56] 


0.12 ±0% 


0.30 ±0% 


0.39 ±0% 


0.20 ±0%f^ 


hetg-11 


4.9 ±8% 


3.5 ±14% 


1.9 ±21% 


[0.54] 


7.9 ±18% 


3.1 ±19% 


[0.56] 


[0.12] 


0.32 ±9% 


0.40 ±5% 


0.20 ±6% 



Note. — Parameter values are followed by an indication of their ± I-ct statistical confidence range, cio to Chi, given as a percentage: 100(^ Chi/cio — 1). A bracketed 
parameter value indicates that it has been fixed at that value during fitting; these fixed values are generally based on other datasets' results. A^i.is is not used (set to 0) for 
the three earliest data sets, due to their lower total source counts. fcTio is fixed for the HETG-99, RGS-03, and the "hetg" datasets which have reduced low-energy counts 
or coverage. fcT^i is fixed for the RGS and HETG-99 datasets because of reduced high-energy coverage and counts. The N abundance is only fit for the RGS and ACIS-00 
datasets, and the O abundance is fixed for the "hetg" datasets which have reduced low-energy coverage and counts. 

''The Nx values given in the table are 10'^ times the XSPEC vpshock models' "norm" parameters, with units of cm~^. Specifically, we have: N^. = 
1x10""'^"'^ (/ n^nHdV) / {iTrD"^) with number densities (rie, nu), volume {dV), and source distance {D) in cgs units. 

''Abundances are relative to the solar photospheric values of I Anders fc Grevessj lll989l V Abundances of Mg and S (not tabulated) are similar to those of Ne and Si, 
respectively. N[j has been fixed at 1.3x10^' cm~^. 
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2.4. The Very-Broad Component and Its Evolution 

To look for and quantify a v-b line component in the SN 1987A spectra we start with 
the data and 3-shock model fits (above) and extend the model to include a v-b component 
that is a broadened version of the 3-shock model itself; the model definition is explicitly 
given in Appendix A, Expression IA2[ The v-b extension adds two additional parameters: 
the fraction of the total flux that is in the v-b component (/v-b) and the width of the v-b 
component expressed as an equivalent FWHM in velocity (fv-b)- 

The choice of broadening the whole spectrum is the simplest given the low signal-to-noise 
ratio of the v-b component in any given line, provided there is some expectation that the 
bright lines of the v-b component are similar to those of the narrow emission. At face value 
this may seem to be unlikely, with the v-b lines produced by a plasma with very different 
temperatures and ionization ages from the plasma producing the narrow lines. However, our 
hydrodynamic modeling shows that the approximation is a surprisingly good one as seen 
in Figure [T3] and discussed in § 14.51 However, based on this comparison we have chosen to 
ignore the FeXVII 0.70-0.75 keV range when doing v-b fitting (below) as it is the only bright 
line-complex that is not predicted to have a v-b component. 

The v-b-3-shock model is fit to the data using ~ 60% finer binning and a somewhat 
smaller energy range that focusses on the brightest, high-resolution lines in the particular 
spectra, generally those of O, Ne, Mg, Si, and Fe-L. For HETG, LETG and RGS data these 
ranges are, respectively: 0.62-2.1 keV, 0.6-1.6 keV (no Si), and 0.47-1.1 keV (including N 
but not Mg and Si). For "hetg-YY" data sets because of their limited low-energy range, the 
O line is not included and the O abundance is frozen at the average of the HETG-07 and 
HETG-11 O v-b values. Because of the reduced energy range of the v-b fitting, the kTi^, 
Tmid and kThi parameters were fixed at their 3-shock values; this leaves the 3 normalizations 
and the abundances to be fit along with the v-b fraction and width. 

As examples, the v-b fits to the HETG-07 and HTEG-11 data are shown in Figure [31 
note that the y-axes of the lower, difference plots are the same, showing that the v-b flux 
is increasing in absolute terms. The statistics of this v-b fitting are given in the left-half of 
Table [3l in particular the values of F2 are large enough to indicate that the addition of the 
two v-b parameters has improved the model fit significantly. 



- 14 - 




Energy [keV] 
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Fig. 3. — Deep HETG spectra and their very-broad (v-b) models. The HETG-07 (top 
pair) and HETG-11 (bottom pair) data are shown (black) with the best-fit v-b 3-shock 
model (red). The non-broad model component is shown as well (blue). The lower half plots 
show the difference (black) between the data and the non-broad component with the v-b 
component over plotted (red). 
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Fig. 4. — Confidence contours for tlie very-broad component parameters. Solid curves give 
tfie 68 % (Ax^ = 2.3) contours for fits to eacfi grating data set; tlie best-fit values are 
indicated with "+"'s. The two deep HETG observations, HETG-07 & HETG-11, show a 
similar, well-determined v-b width of ~9300 kms~^ FWHM. See the text for a discussion of 
the spread seen among the contours. 



Instead of tabulating the best-fit v-b parameters and their ranges, it is more useful to 
generate two-dimensional confidence contours in "v-b fraction vs. v-b width" space. These 
contours and their best-fit values are shown in Figure HI for HE TG-99 the dat a poin t with 



errors is based on the 2-Gaussian "stacked fitting" result given in lMichael et al.l ( 12002| ). The 
earliest data sets, HETG-99 and RGS-03, show large values of v-b fraction (> 0.6) with 
widths of order 10x10^ kms^^. In contrast the later observations show v-b fractions within 
the range 0.15 to 0.30 and cover a large spread in widths, with the HETG values below 
15x10^ kms~^ and RGS-determined widths above this value. 

This separation suggests that the measurement of the v-b component depends on the 
particular spectrometer used. Figure E] shows how the spectrometers "see" a monochromatic 
line consisting of a narrow and a v-b-component; these plots include SN 1987A's spatial- 
spectral blur descibed in Appendix |X1 The HETG, even at the higher MgXII energy, clearly 
separates the narrow and v-b components, although it's also clear that the v-b presence is 
subtle. For the LETG the narrow component is significantly broadened by the spatial extent 
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of SN 1987A and at high energies the v-b component is almost completely covered by the 
narrow component. The RGS line-spread function is more peaked than the LETG at OVIII, 
however it includes extended wings so that even the narrow component produces counts far 
from the line center. Hence the RGS v-b measurement will be particularly sensitive to the 
accuracy of its LSF calibration. 

With these considerations in mind, we see that the two deep HETG observations, HETG- 
07 and HETG-11, should give the most accurate measure of the v-b component. Comparing 
their small confidence contours, there is a clear change in the v-b fraction from 2007 to 
2011. In terms of the v-b width, these two HETG contours show very similar values of ~ 
9300 kms"^ FWHM, with an estimated 90 % confidence range of ±2000 kms ^. This value 
is also in rough agreement with the (lower statistics) HETG-99 and RGS-03 widths and so 
it is reasonable to postulate a constant v-b width with a fraction that changes over time as 
shown by the blue path in Figure IH 

Given the above, we fix the width of the v-b component at v^-h = 9300 kms~^ FWHM 
and fit for just the v-b fraction. The results of this "constrained" fitting are given in the 
right-half of Table 121 including the v-b fractions and their 68 % confidence ranges. Note 
that for most data sets the change when doing this constrained fitting, Axf, is a sub- 
stantial portion of the change when the v-b width is also being fit, Ax2- This suggests 
that the data do not strongly constrain the width to differ from our assigned value; the 
elongated contours of Figure H] suggest this as well. Conversely, the v-b fractions given by 
the constrained fitting have only a small dependence on the assigned width. As examples, 
when the data sets LETG-07, RGS-11, HETG-11, and hetg-11 were fit with the width set to 
6500 & 13000 kms-i FWHM, the best-fit fractions found are: 0.132 & 0.160, 0.216 & 0.210, 
0.163 & 0.178, and 0.200 & 0.188, respectively. These variations are actually within the sta- 
tistical 1-0" ranges given for the /v_b values in Table |31 and so the choice of the v-b width 
will not significantly change the v-b fractions and light curve. 

Having determined the v-b fractions for the data sets we can now re-do the 3-shock 
fitting of § 12.31 with the appropriate v-b component fixed and included in the model. This 
gives the results in Table |H Note that the abundances have increased compared to the 
non-v-b values in Table |21 This is expected since a larger flux in a given line is now needed 
in order to match the height of the narrow-line component in the high-resolution spectrum. 
The change is especially pronounced for HETG-99 where the v-b fraction is 0.78 and the 
abundances are now in better agreement with those of ACIS-00 (although there is a large 
error range.) Also as expected, the ACIS-00 values hardly change when the v-b component 
is included because the lower-resolution CCD spectra are primarily determined by the total 
line flux with little sensitivity to any underlying velocity structure. 
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Fig. 5. — Comparing the grating spectrometer line responses when observing SN 1987A. 
The HETG(MEG), LETG, and RGS views of the same combination of a narrow hne plus a 
very-broad component are shown for the case of /v-b = 25% and Wv-b = 9300 kms~^ FWHM. 
Because the resolution varies with energy, responses are shown for O VIII (left colum) and 
for Mg XII (right column) . On each plot the orange curve shows a continuum level, the blue 
curve is the continuum plus the v-b component, the magenta curve is the continuum plus 
the narrow component, and the red curve is the sum of all three. Vertical gray lines indicate 
the line center and ±1500 kms~^ from it. For reference, example data are shown in the 
background by the light gray histograms; the OVIII region includes the He-like L-/3 line at 
~ 0.665 keV. The HETG provides the best separation of the v-b and narrow components 
especially at the higher energies. 
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Table 3. Very-Broad Fitting Statistics and Results 



Data ID 






b 


Aa 


■2 


F2 


Ax 


,2 


-b = 9300 kms- 

i^l /v-b ■ 


-1 FWHM — 
l-o" range 


HETG-99 " 


32- 


9 


24, 


,0 


2, 


,0 


1, 


,0 


1. 


,2 


1. 


,1 


0. 


,372 


0, 


,02 - 


- 0. 


,63 


RGS-03 


226- 


-8 


352, 


,7 


86, 


,0 


26, 


,6 


80. 


,0 


48. 


,9 


0, 


,655 


0, 


,61 - 


- 0. 


,70 


LETG-04 


130- 


-9 


174, 


,3 


68, 


,6 


23, 


,8 


66. 


,9 


46. 


,4 


0, 


,248 


0, 


,22 - 


- 0. 


,28 


RGS-07 


216- 


-9 


360, 


,5 


36, 


,3 


10, 


,4 


31. 


,0 


17. 


,7 


0, 


,175 


0, 


14 - 


- 0. 


.20 


HETG-07 


622- 


10 


1146. 


1 


347, 


,8 


92. 


,9 










0, 


,244 


0, 


,23 - 


- 0. 


,26 


LETG-07 


140- 


-9 


355, 


,2 


81, 


,7 


15, 


1 


48. 


,6 


16. 


,5 


0, 


,150 


0, 


,13 - 


- 0. 


.17 


RGS-08 


231- 


-9 


450, 


,3 


71, 


,6 


17. 


,6 


45. 


,7 


21. 


,4 


0, 


,181 


0, 


,16 - 


- 0. 


,21 


hetg-08 


146- 


-9 


174, 


,4 


28, 


,9 


11. 


,4 


26. 


,4 


20. 


,6 


0, 


,254 


0, 


,21 - 


- 0. 


,30 


RGS-09 


217- 


-9 


416, 


,1 


68, 


,9 


17. 


,2 


33. 


,7 


15. 


,6 


0, 


,158 


0, 


,13 - 


- 0. 


.18 


hetg-09 


349- 


-9 


546, 


,6 


76, 


,3 


23. 


,7 


75. 


,0 


46. 


,7 


0, 


,208 


0, 


.18 - 


- 0. 


.23 


RGS-10 


222- 


-9 


404, 


,0 


81, 


,2 


21. 


,4 


66. 


,8 


34. 


,2 


0, 


,206 


0, 


,18 - 


- 0. 


,23 


hetg-10 


328- 


-9 


427, 


,0 


44, 


,5 


16. 


,6 


42. 


,4 


31. 


,6 


0, 


,178 


0, 


,15 - 


- 0. 


.20 


RGS-11 


201- 


-9 


429, 


,6 


67, 


,1 


15. 


,0 


57. 


,5 


25. 


,2 


0, 


,211 


0, 


,19 - 


- 0. 


.24 


HETG-11 


599- 


10 


1115, 


,3 


116, 


,7 


30. 


,8 










0, 


,173 


0, 


16 - 


- 0. 


19 


hetg-11 


319- 


-9 


422, 


,0 


47, 


,2 


17. 


,3 


45. 


,1 


33. 


,1 


0, 


,201 


0, 


,17 - 


- 0. 


,23 



Note. — The four left-most numerical columns give values related to the very-broad 
fitting of the data sets: the number of data bins minus the number of free fit parameters, 
the value for the v-b fit, the amount by which Xv-b improvement over (less than) 

the non-v-b value, and the F-statistic for adding the 2 v-b degrees of freedom to the 
model. Specifically, F2 = (Ax^/S) / (Xv-b/ (A^data ~ ^free)) a-nd values of F2 greater than 
_F^j.it ^ 7.3 indicate a high confidence, > 0.999, that the addition of the 2 v-b parameters 
does result in a significantly improved model. 

The right-most four columns give values related to a constrained v-b fitting where the 
v-b velocity width is fixed at the HETG-07/HETG-11 average value; thus there is only a 
single free v-b parameter: /v-bi the v-b fraction. The first two of these columns give the x^ 
improvement and the appropriate i^-statistic: Fi = (Axf /l)/(x93oo/(-^data ~ Afj^^ge -|- 1)) 
with X93Q0 ~ Xv-b Axi ~ Axi- In this case we have F^rit ~ 11.4. The right-most two 
columns give the best-fit v-b fraction and its l-cr confidence range for the constrained 
fitting. 

^The HETG-99 values tabulated here are from the same analysis as us e d for the other 
data sets in this table. The "stacked analysis" result of iMichael et all |20o2)> /v-b = 
0.78 ± 0.1, is considered more accurate and is used for the very-broad HETG-99 point in 
Fig.[TD] 



Table 4. Parameters of the Global 3-Shock-witli- Very-Broad-Component Model Fits 







Normalizations 


a 


fcTio 


mid 






— Relative abundances ^ — 




Data ID 








(keV) 


(10"^) 


(keV) 


N 


O 


Ne 


Si 


Fe 


HETG-99 


0.10 ±54% 


[0] 


0.12 ±9% 


[0.74] 


9.2 ±211% 


[4.2] 


[0.56] 


0.82 ±214% 


0.77 ±154% 


0.45 ±120% 


0.26 ±144% 


ACIS-00<= 


0.18 ±40% 


[0] 


0.12 ±9% 


0.73 ±13% 


8.2 ±69% 


4.2 ±17% 


1.6 ±142% 


0.47 ±25% 


0.61 ±73% 


0.68 ±43% 


0.12 ±24% 


RGS-03 


1.5 ±21% 


[0] 


0.17 ±91% 


0.57 ±8% 


3.0 ±24% 


[4.3] 


0.86 ±29% 


0.16 ±23% 


0.28 ±16% 


0.68 ±145% 


0.11 ±25% 


RGS-07 


2.4 ±7% 


0.16 ±64% 


1.0 ±10% 


0.52 ±2% 


5.2 ±8% 


[3] 


1.6 ±7% 


0.21 ±6% 


0.64 ±8% 


0.81 ±15% 


0.36 ±6% 


RGS-08 


2.6 ±2% 


0.40 ±21% 


1.2 ±6% 


0.53 ±2% 


5.7 ±1% 


[2.7] 


1.9 ±3% 


0.25 ±4% 


0.62 ±7% 


0.54 ±20% 


0.43 ±4% 


RGS-09 


4.3 ±3% 


0.86 ±9% 


1.4 ±9% 


0.54 ±1% 


8.1 ±8% 


[2.7] 


1.6 ±5% 


0.23 ±2% 


0.57 ±6% 


0.44 ±18% 


0.33 ±4% 


RGS-10 


3.8 ±4% 


0.98 ±15% 


1.6 ±9% 


0.55 ±1% 


12 ±3% 


[2.7] 


2.5 ±4% 


0.33 ±3% 


0.80 ±5% 


0.81 ±10% 


0.42 ±3% 


RGS-11 


3.4 ±2% 


0.68 ±16% 


2.1 ±8% 


0.57 ±2% 


17 ±8% 


[2.7] 


3.3 ±2% 


0.49 ±1% 


1.2 ±5% 


0.93 ±14% 


0.59 ±0% 


LETG-04 


1.8 ±6% 


0.41 ±11% 


0.34 ±7% 


0.50 ±3% 


4.3 ±4% 


4.5 ±14% 


[0.56] 


0.14 ±10% 


0.38 ±7% 


0.40 ±7% 


0.18 ±5% 


HETG-07 


3.9 ±2% 


0.60 ±3% 


1.1 ±1% 


0.50 ±1% 


3.4 ±2% 


3.0 ±1% 


[0.56] 


0.15 ±4% 


0.33 ±2% 


0.41 ±2% 


0.22 ±2% 


LETG-07 


4.6 ±2% 


0.84 ±3% 


1.5 ±4% 


0.49 ±2% 


3.8 ±2% 


2.3 ±3% 


[0.56] 


0.11 ±3% 


0.33 ±2% 


0.39 ±3% 


0.22 ±2% 


hetg-08 


5.6 ±9% 


0.54 ±189% 


1.9 ±25% 


[0.52] 


4.1 ±15% 


2.2 ±15% 


[0.56] 


[0.15] 


0.32 ±9% 


0.42 ±6% 


0.19 ±11% 


hetg-09 


4.5 ±6% 


1.2 ±23% 


1.9 ±10% 


[0.52] 


5.5 ±13% 


2.6 ±6% 


[0.56] 


[0.15] 


0.41 ±8% 


0.46 ±6% 


0.28 ±6% 


hetg-10 


5.2 ±7% 


1.9 ±22% 


2.2 ±15% 


[0.52] 


5.7 ±17% 


2.6 ±9% 


[0.56] 


[0.15] 


0.35 ±9% 


0.46 ±6% 


0.23 ±7% 


HETG-11 


5.3 ±2% 


2.2 ±1% 


2.2 ±1% 


0.55 ±1% 


7.2 ±2% 


2.8 ±2% 


[0.56] 


0.15 ±7% 


0.40 ±2% 


0.49 ±2% 


0.26 ±2% £ 


hetg-11 


4.1 ±9% 


2.4 ±19% 


2.3 ±15% 


[0.55] 


7.7 ±29% 


2.8 ±11% 


[0.56] 


[0.15] 


0.45 ±14% 


0.51 ±9% 


0.26 ±7% 



Note. — The data sets have been fit as before, Table[2l with the inclusion of a v-b component fixed to have FWHM = 9300 kms^^ and a fraction as given in TablefS] 
Parameter values are followed by an indication of their ± l-cr statistical confidence range, cio to Chi, given as a percentage: 100(y^cj^i7cio ~ !)• A bracketed parameter 
value indicates that it has been fixed at that value during fitting; these fixed values are generally based on other datasets' results. A^i.is is not used (set to 0) for the three 
earliest data sets, due to their lower total source counts. kTi^ is fixed for the HETG-99 and the "hetg" datasets which have reduced low-energy counts or coverage. fcTjj; is 
fixed for the RGS and HETG-99 datasets because of reduced high-energy coverage and counts. The N abundance is only fit for the RGB and ACIS-00 datasets, and the O 
abundance is fixed for the "hetg" datasets which have reduced low-energy coverage and counts. 

''The Nx values given in the table are 10'^ times the XSPEC vpshock models' "norm" parameters, with units of cm~^. Specifically, we have: Nx = 
1x10^^^ (/ neTLjjdV)/ {'i-irD'^) with number densities (ne, n^f), volume {dV), and source distance (D) in cgs units. 

''Abundances are relative to the solar photospheric values of I Anders fc Grevessj lll989l) . Abundances of Mg and S (not tabulated) are similar to those of Ne and Si, 
respectively. Njj has been fixed at 1.3x10^' cm~^. 

'^A v-b fraction of 0.7 was used in fitting the ACIS-00 data here; the fit results are very similar to those for ACIS-00 in Table [2] (where v-b fraction = 0). 
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HYDRODYNAMIC MODELING OF SN 1987A 



Besides being an extremely well-studied object, SN1987A has been frequently mod- 
elled in almost every wavelength regime. The earliest calculations, which delineated the 
circumst ellar material (CS M) into a low-density inner and high-de nsity outer r e gion, were 
made by lltoh et al.l (119871) . The X-ray emission was modelled by iMasai et al.l ( 119911 ) and 



Luo fc McCravl ( 
were explo red by 



991ar) using 1- d imensi onal (ID) models. Two-dimensional (2D) models 



Luo fc McCravl (Il991bl ) for the case of an hour-glass equatorial ring (ER) 



structure. ISuzuki et al.l (119931 ) carried out 2D axisymmetric smoothed particle hydrody- 
namics (SPH) simulations with the ER modeled as a torus; their Figures 4-8 demonstrate 
the hydrodynami c comp lexity by showing the different components' spatial distributions. 
Masai &: Nomotol (119941 ) went back to the ID models but t ook into account the light-travel 
times from different parts of the remnant. iLuo et al.l (119941 ) carried out 2D calculations ap- 
proximating the ER as a rigid boundary; they also introduced the ejecta density distribution 
that has been subsequently used by most authors, and forms the basis for the ejecta density 
profile used in this paper. 



The increasing X-ray and radio emission led IChevalier fc DwarkadasI (119951 ) to propose 
the existence of an ionized Hll region interior to the ER. The FS had collided with this 
region around day 1200, and was expanding in this moderately-high density medium, p ^ 
100 amucm"^, giving rise to increasing X-ray and radio emission. The slowing down of 
the shock wave observed at radio wavelengths was consistent w i th this conjecture. The X- 
ray flux from the Hll region was modelled by iBorkowski et al.l ( Il997al ). whereas a detailed 
prediction of the expected X-ray fl ux and spectrum when the shock collided with the dense 
ER was given in IBorkowski et al.l (Il997b( ) . 



For the next decade or so, investigations of the X-ray emission were mainly made via 
simpler, non-hydro dyn amical models that attempted to tie in the observed X-ray emissio n 

to the CSM proper ties JPark et al.lbooi boosi bood . Iflaberl et allbood : Izhekov et aPboioh. 

DwarkadasI (l2007al ) revised the older simulations made (but not presented) in I Chevalier fc Dwarkadas 
(|l995l ) using more recent data, and managed to create spherically symmetric hydrodynam- 
ical models that successfully reproduced the FS radius and velocity as measured at radio 
frequencies. Using these simulations he w as also able to make a crude estimate of the X-ray 
light curve. This model was updated in iDwarkadad (l2007bl ). and forms the basis of the 
hydrodynamical simulations that are presented in the following. 
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Fig. 6. — Schematic diagrams of the CSM as modeled in lZhekov et al.l (120101 ). top row, and in 
this work, bottom row. The left diagrams show the CSM configuration assumed at the time of 
the SN explosion. The coresponding diagrams at the right show the configuration ~ 7300 days 
after explosion when the forward shock, driven by the expanding ejecta, has progressed into 
the dense, clumpy ER. The Zhekov model is based on a single ID hydrodynamics whereas 
our model consists of two ID models, one that includes the dense ER and one that accounts 
for the Hll material above and below the equatorial plane. It is this latter component that 
gives rise to the very-broad X-ray emission. 



3.1. Our Hydrodynamic Model of SN 1987A 

The structure of the CSM around SN 1987A is inherently 3-dimensiona l, as is visible 



i n the amazing images from the Hubble Space Telescope, most recently in iLarsson et al. 



(I2OIII ). Modelling the ejecta-CSM interaction in 3-dimensions is a highly challenging and 
computationally demanding problem, does not encourage iterative parameter exploration, 
and requires multi-dimensional non-equilibrium ionization (NEI) X-ray emission calcula- 
tions. However, given the substantial symmetries in the overall geometry, it is possible to 
partition the domain into one or more regions whose emission can then be approximated by 
a spherically symmetric ID simulation. These simulations can then be summed together, 
weighted by the fraction of the full in solid angle that they each subtend. 

As examples o f this approach we show schematics in Figure [6] for the modeling of 



Zhekov et al.l (120101 ) and our model. The former consists of a single ID solution with a 



smooth, increasing-density CSM component. To match the observed spectra, especially at 
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late times, dense "clumps" are included in the CSM with the c lump X-ray emission based on 



transmitted and reflected shock equations (jZhekov et al.ll2009l ). Our model in the lower por- 



tion of Figure E] is motivated by the on-going existence of the v-b component which requires 
a combination of two ID solutions: 

i) A "with ER" simulation that consists of an Hll region followed at larger radii by a 
density jump to high values at the ER location; the density proflle is shown in Figure [71 
This component subtends a small solid angle, ^Ier, and the high-density ER includes 
additional "c l umpin g" (§ I4.2l fc l47i|) . This component is similar to the ID model of 



Zhekov et al.l (|2010|) 



ii) A "no ER" simulation in which the Hll region continues to large radii and there is 
no ER density jump. This simulation produces the v-b emission and corresponds to 
material above and below the equatorial plane; it subtends a moderate solid angle, 

^HII- 

This multi-ID technique does have its limitations, especially for the "with ER" case, and we 
will keep these in mind. 

The simulations here were carried out using the VH-1 code, a 1, 2, and 3-dimensional 



finite -difference hydrodynamic code based on the Piecewise Parabolic Method (jColella fc Woodward 



19841 ). The code solves the hydrodynamical equations in a Lagrangian framework, followed 
by a remap to the original Eulerian grid. Radiative cooling is included in the form of the X- 
ray cooling function scaled by a factor of 3.5 to account for IR dust emission, Ljr/Lx ~ 2.5. 

For each of a set of time steps, the code outputs a snapshot file giving the hydrodynamic 
quantities density, velocity, and pressure at each radial zone; two examples are shown in 
Figure [HI The set of snapshot files is then processed by custom routines coded in S-Lang0 
and run within the ISIS software environment. The radial grid is rebinned into mass cuts, 
so that the time evolution of each parcel of gas in a given mass shell can be accurately 
followed; in this way t he NEI state is tracked and the X-ray emission can be calculated 



(iDwarkadas et al.ll2010l ). 



Of course we did not simply pick parameters, run the simulation, and have agreement 
with observations. An iterative process of trial & error was carried out to converge on model 
parameters that reproduced the various measured properties of SN 1987A, especially in the 
X-ray. Rather than take the reader through these steps (and mis-steps), we summarize in 



^ S-Lang web page: http : //www, jedsof t . org/slaingTI 
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the following sections the final model properties and then compare the model results to data 
in §11 



3.2. Ejecta Parameters 



The ejecta are modelled using the IChevalierl ( 1l982l ) prescription where the ejecta are in 



homologous expansion with a density that decreases as a power-law with radius, i.e., 

Pej = C t-^ {r/ty- . (1) 

This is written to show that one can also consider the density profile as a function oi v = r /t. 
Since a power-law extending all the way back to the origin is unphysical, below a certain 
velocity, Vt, the density is assumed to be a constant, see the "Ejecta" portion of the density 
in Figure [71 Thus the profile can be specified by the three parameters: C, n, and Vt . The 



eject a's total kinetic energy and mass are functions of these parameters (jChevalier fc Liang 



19891 . equ. 2.1) and so we can equivalently choose other 3-parameter sets, e.g.: n, i^ej, and 



In our ID hydrodynamic simulations, the RS moves into the ejecta and at most shocks 
only the outer ~ 0.5 solar masses of the (47r) ejecta profile, at which point the ejecta velocity 
is ~5700 kms~^. Provided that Vt is less than this value, which is generally the case, we 
are insensitive to the plateau transition location. Our simulations then only depend on the 
outer profile, that is the two parameters of Equation [T] the norm alization of the power-la w 



profile (C =2.03x10^^ gcmS-6)l3 and the exponent = 9 from ILuo fc McCravl fll991bf )l 
Implicitly fixing these two, we are then left with a degenerate degree of freedom, our insen- 
sitivity to ff, in choosing the ejecta parameters. Hence there are pairs of (£'ej,Mej) that will 
equivalently give the hydrodynamics presented here: (2.4,16), (1.5,7.9), or (1.0,4.3), where 
the energy is in units of 10^^ erg and the mass is in solar masses. It's important to note that 
we do not make any determination of the actual global values of these parameters, they are 
only used to define the outer ejecta properties relevant to our hydrodynamics. 



^ The value and units of C do not mean much at face value; following many authors we can instead specify 
the value of poj at some appropriate reference values of v and t. For w = Ix 10^ kms~^ and t =10 y ears we 



get pcj ~ 390 amu cm ^ , very similar to the value of 360 amu cm ^ used by iBorkowski et al.l (|l997al ) 
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3.3. Hll Region Parameters 



Estimates of the paramet ers of the Hll re gion are determined mainly from observa- 
tions. This was first done by iLundqvistI (Il999l ) by analyzing the ultraviolet line emission 
from SN 1987A. iDwarkadad (l2007al Jbl) attempted to refine these parameters by calculating 
the radius and velocity of the expanding SN shock wave in a spherically symmetric case, 
comparing the results to radio observations, and iterating until a good fit was obtained. 
They also calculated a reasonable fit to the hard X-ray emission under collisional-ionization 
equilibrium (CIE) conditions. 

In the present work we further develop these earlier calculations. The H II region radius 
and density profile are adjusted so that our simulated emission agrees with the early X-ray 
light curves (up to day ~5000); in particular, producing a reasonably sharp "turn-on" at 
~ 1400 days and agreeing with the HETG-99 and ACIS-00 measured spectra at later times. 
A further constraint on the H II region, especially near the equatorial plane, comes from the 
requirement that the FS reaches the ER (an d/or protrusions) at a time and radius in agree- 
ment with the X-ray imaging observations (IRacusin et al.l 120091 ). The parameters we have 
found appropriate for the Hll region, shown in Figures [71 have it beginning at 3.61x10^^ cm 
(0.117 pc) with a density of ~ 130 amucm"'^ and gradually increasing to ~250 amucm^^ 
at 6.17x10^^ cm (0.20 pc) and remaining at this density to larger radii. These X-ray-based 
values have an inner radius closer to t he expectations of I Chevalier fc Dwarkadad (119951 ) 
and differ from the previous Hll values (iDwarkadasI l2007bl . beginning at 4.3x10^^ cm with 
p ^ 6 amucm"^, increasing to p ~ 200 amucm"^ at 6.2x10^^ cm) which were based primar- 
illy on the radio size at early time; we discuss the reconciliation of these further in § 14.11 and 

§S21 



3.4. Equatorial Ring Parameters 



Unlike the Hll region where the simple ID approximation may suffice, it is clear that 
only a far cruder approximation to the X-ray emitting ER can be made with simple ID 
constructs. Note that in this paper we often use the term "ER" as a generic term to refer to 
the dense (~ 10^ amucm"^) equatorial material with which the SN shock interacts. The ER 
therefore includes the finger-like extensions that iMcCrayl (120071 ) has postulated as extending 
inwards from the ring. Consequently, the ER is not located at any sin gle radius; this is clearl y 
demonstrated by the gradually increasing number of optical hot spots (ISugerman et al.ll2002l ). 
A superposition of ER collisions in time is thus expected to make up the ER contribution 
to the X-ray light curve. 
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Radius (pc) 

Fig. 7. — Initial ID model radial profiles. This early-time (~ 900 d) configuration shows 
the ejecta, the Hll region, and the options of a "thick" or "thin" ER. The ejecta (black) 
and CSM (blue) densities are separated at the contact discontinuity (CD); the mean plasma 
temperature (red) and velocity (green) are also plotted. The gray vertical hues indicate the 
boundaries of the mass-shells that are tracked and used to compute the NEI X-ray emission; 
note that they do not need to extend into the plateau region of the ejecta profile. 



For modeling purposes we choose a single location of the ER starting at a radius of 
5.4x10^^ cm (with a density of ~9000 amucm"'^). As we'll see, this choice reproduces the 
"kink" in the measured X-ray radius vs time (§ 14. ip and the bulk of the dramatic soft X-ray 
increase after day 6000 (§ 14. 2p : hence this model component is representative of the majority 
of shocked ER na aterial. Since there is an indication that the X-ray flux may be leveling off 
( jPark et al.ll201ll ) we also consider the case of a finite or "thin" ER with a density that then 
drops from ~ 10800 amucm"^ to (a lower, arbitrary value of) ~ 1200 amucm"^ at a radius 
of 5.77x10^^ cm, as shown in Figured 
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Fig. 8. — Model profiles at later times. Top: The simulation at the ACIS-00 epoch when the 
FS has moved well into the Hll region; note that the early-time density profile is shown in 
gray for reference. Bottom: At the HETG-11 epoch, the FS has encountered the "thin ER", 
slowed, and is just about to exit the ER. The color coding is the same as for Figure [71 



4. MODEL RESULTS AND COMPARISONS 



The previous section summarized the ejecta-CSM properties of our hydrodynamic sim- 
ulations. In the following we compare several results derived from the simulations with the 
observed properties of SN 1987A. 



4.1. Model Locations 



The radio and X-ray measured radii are used to constrain the hydrodynamic models. In 
Figure [9] we show these radii along with the locations of the simulations' FS, CD, and RS. 
The measured r adii are given i n arc second s and are based on d e-projected model fitting for 
both the radio ( iNg et al.ll2008l ) and X-ray (iRacusin et al.ll2009l ) values. To convert from the 
simulation spatial units (cm) to observed angular units a distance of 48.5 kpc has been used 
to get the agreement shown; this is very close to the accepted distance to SN 1987A (§[1]) 
and further "fine tuning" of the hydrodynamics is not warranted for these simple models. 

For the "with ER" simulation, Figure M top panel, the locations of the model emitting 
region (between the FS and CD) are in good agreement with the measured ACIS radii 
and reproduce the " kink" that occurs at about day 6000 at a radius of ~0.74 arc seconds 
(IRacusin et al.ll2009[ ). This radius is set by the start of our ER at 5.4x10^^ cm, or 0.744 
arc seconds, and it also agrees well with the inner radii of the majo r ity of the optical ring 
material that is seen in the very useful Figure 7(b) of ISugerman et al.l (|2002| ). Hence, our ER 
model location is a good representation of the location of the start of the bulk of the ER. 
Regarding a possible end to the ER, our "thin ER" case ends at a radius of 5.77x 10^'' cm, or 
0.795 arc seconds. This gives our thin ER a radial extent of ~ 0.05 arc seconds, and, looking 
again at Sugerman's Figure 7(b), we see that this is a reasonable radial extent for the central 
portions of individual bright optical features. At late times, the locations shown in Figure M 
are based on the "thin ER" profile (labeled in Figure[7]) which has the FS leaving the dense 
ER at ~ 9000 days; this produces a slight upward kink in the FS location curve at this time. 

The "no ER" simulation, representing the out-of-plane Hll region, does a good job of 
matching the ACIS radii before day 6000 (since it is the same as the "with ER" simulation 
at those times) and continues to expand at a roughly constant rate. By 25 years, the 
FS in the Hll region simulation has progressed to a distance of ~ 6.8x10^^ cm (0.22 pc), 
equivalent to ~0.94 arc seconds for the model-scaling distance of 48.5 kpc. Comparing this 
to the mean radius and extent of the optical ER, 0.83 and 0.65-1.0 arc seconds, respectively 
(ISugerman et al.ll2002l . Figure 7), suggests that the FS in the out-of-plane Hll region is now 
beyond most of the ER. 
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Comparing the model locations with radio measurements on these plots, the 8.6 GHz 
data (INg et al.ll2008l ) appear to be different from both the "with ER" and "no ER" X-ray- 
based models; in particular showing a > 0-6 ar c second radius by day 2000. However, using 
the previous CSM profile of iDwarkadad (l2007bl ) does give a CD-FS location curve that is in 
good agreement with the measured average radii for the 8.6 GHz radiat ion, Figure |9] lower- 
right panel. The single 36 GHz radius measurement ( jPotter et al.l 120091 ) is located between 
the X-ray and the 8.6 GHz values, although it has an uncertainty range including each of 
these. One way to accomodate the larger 8.6 GHz radio radii is to have that emission come 
from material further from the equatorial plane as schematically shown in Figure [151 



4.2. Model Light Curves 



We use the hydrodynamic simulations as input to calculate the X-ray emission from 
each shocked mass-cut shell at each time step, employing the s ame scheme to tr a ck th e 
NEI plasma state and carry out spectral calculations as outlined in lDwarkadas et al.l ( 120101 ). 
Once the hydrodynamic solution is determined - fixing temperatures, ionization ages, and 
normalizations for the mass-cuts - the only free parameters are the abundances, an optional 
clumping set by a clumping fraction and its over-density, and the overall normalization set 
by the solid angle which each ID model actually fills, flnu and ^er. The abundances 
can be determined by fitting to the measured spectra, as described further in § l4.3l fc l4^ 
In terms of line broadening, the hydrodynamic velocities of the shells are used to include 
Doppler broadening in the synthesized spectra; in § 14.51 this "real" broadening is compared 
with the simple smoothing approximation that was used to quantify the v-b fraction (§ 12. 4p . 
Finally, for high shock velocities and low r's, the simple oc f ~^ behavior oi = Tp/T ,; 



s in 


Dwarkadas et al. 


(2010) 


van Adelsbere; et al. 


2008^ 



to set a minimum value of /?„ 



0.06. 



Integrating the simulated spectra over the standard energy bands, 0.5-2 keV and 3- 
10 keV, we created model light curves, shown in Figure [TOl The light curves for the two 
hydrodynamic simulations have been plotted separately and their (imagined) sum can rea- 
sonably describe the observed SN 1987A light curves over the full time span. 
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Fig. 9. — Model radii compared with X-ray and radio measurements. Top: The "with ER" 
simulation has most X-ray emission between the FS (blue) and CD (red) radii and their 
locations agree very well with the ACIS-measured average radii (black). Note how the kink 
due to the ER collision occurs first in the FS curve (around day 6000) and then at later times 
in the CD and RS (green) curves. Lower-left: As expected, the "no ER" simulation shows 
continued expansion; this applies to the case of the out-of-plane H II region giving rise to the 
very-broad X-ray component. Lower-right: The 8.6 GHz average radii (orange) are b etter 
matched with a difi^erent CSM profile, here we use the one given in iDwarkadad (l2007bl ). 
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Days since SN 1987A 



Fig. 10. — Comparing X-ray light curves with the hydrodynamic model fluxes. The Hll 
region emission from the "no ER" model is shown by the curves beginning before day 3000 
(magenta, green). The model does a reasonable job of matching the initial and late-time 3- 
10 keV flux (green diamonds) as well as modeling the initial 0.5-2 keV turn-on and growth 
(ROSAT points). At late-times the 0.5-2 keV flux from the Hll region is in reasonable 
agreement with the measured very-broad flux (magenta diamonds), strengthening the iden- 
tification of the v-b emission with the Hll region. The emission from the shock collision 
with the dense ER rises steeply around day 6000 (orange, light blue). The light curve here 
(orange) is for the the "thin (or finite) ER" case of Figure[71 the gray line indicates the 
0.5-2 keV flux when the ER continues in the "thick ER" case. 
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A filling fraction of Qhu/Att = 0.26 is used to match the overall flux level for the 
H II ( "no ER" ) simulation, corresponding to emission from ± 15 degrees of the ring plane. 
This is roughly simi lar to other esti mates of the H II reg ion's out-of-plane e xtent, e.g., by 



Michael et all fcOOl ±30 degrees), IZhekov et al.l fcOlOl . ±10 degrees), and iMattila et al. 



J2OI0I . ±7.2 degrees). The abundances for the Hll light curve are set from fits to the early 



ACIS observation, described in § 14.31 One characteristic of this early-time H II emission is its 
relative hardness, e.g., for the deep ACIS-00 observation at day 5036 the ratio -F3-io/-^o.5-2 
is 0.35 whereas at later times (day 7997, January 2009) the ratio has dropped to 0.12. As 
the light curves show, this is caused by the increase in the 0.5-2 keV flux due to the collision 
with the ER, rather than any decrease in the 3-10 keV emission. 

In contrast to the Hll region, the ER needs to produce emission primarily in the 0.5- 
2 keV range while keeping the 3-10 keV flux below the measured levels. Even though the ER 
is dense, our simulation still gives electron temperatures at/above 3 keV and thus relatively 
too much 3-10 keV flux. We could further increase the density of the ER, however given the 
resonable agreement of the slope of the "with ER" simulation (top panel of Figure |9]) and the 
expectations of an inhomogeneous ER (e.g., the optical spots), we consider another, approx- 
imate approach. We can match the observed spectra by including emission from "clumps" 
in the ER with a density enhancement ^^ciump ~ x5.5 that fills /dump ~ 30% of the CSM 
volume. The clumping is implemented at the spectral synthesis stage and introduces kTe ~ 
0.5-0.8 keV emission components greatly enhancing the 0.5-2 keV emission (further details 
of the clumping implementation are given in § 14.41 ) We acknowledge that for these clumping 
parameters most of the mass is in the very dense clumps (i.e., 0.30x5.5 / (0.70x1+0.30x5.5) 
~ 70%) and so the clumping is not a small perturbation on the ER hydrodynamics; there- 
fore our ID clumped ER simulation is just a placeholder for a realistic multi-dimensional 
simulation. 

The single set of "with ER" light curves plotted in Figure (TU] represents the bulk of the 
0.5-2 keV flux increase; the fractional area used to scale the ID model's emission is very 
small: flER/4:7i = 0.0016 . To put this into perspective, the fraction subtended by a uniform 
region with the dimensions of the optical ring, i.e., a region extending ±0.05 arc seconds 
out of the plane at a radius of 0.83 arc seconds, is almost 40x larger, ~ 0.0600. Hence, the 
small value indicates that the the ring is not uniformly dense and clumped, but rather it has 
"X-ray hot spots" analogous to the (still) discrete optical spots around the ring. Specifically, 
the 0.0016 value would arise from ~20 regions each subtending a diameter of ~2 degrees 
around the ring. 

It is seen in Figure [10] that there is earlier-time flux growth above the H II emission 
that our plotted ER curve does not include. As mentioned in § I3.4[ a range of inner radii is 




Fig. 11. — Left: The distribution of kT^ and r values in the Hll hydrodynamics at the ACIS- 
00 epoch. Points of (r, kTe) are plotted for all shocked mass shells in the model; symbols for 
the shocked ejecta (black) and CSM (blue) are larger/bold for brighter shells. Right The 
ACIS-00 data (black) are fit with a model (red) that is the sum of the emission from the 14 
shocked shells of the left panel; a single norm and the abundances of N, O, Ne, Mg, Si, S, 
and Fe were the only free parameters in the fit. The emission from the CSM shells (blue) is 
much larger than that of the shocked ejecta (green) at this early epoch. 



expected for the "ER" with some material at radii smaller than the 0.74 arc seconds where 
our modeled ER starts. Mentally "translating" our ER curve, we would expect to get flux 
beginning ~5500 days by having a small fraction of ER material (perhaps 0.1 x Qer) that 
is impacted at ~4600 days. From our hydrodynamics this time corresponds to a radius of 
4.8 xlO^^cm (0.155 pc) or 0.67 arc seconds. This radius is in reasonable agreement with the 
locations of the early optical spots (jSugerman et al.ll2002l . Figure 7). 



4.3. Hll Abundances from the Hydrodynamics 



Using the set of shocked mass-cut shells at a given time-step of the hydrodynamics, we 
can create a spectral model in ISIS which is the sum of emission from each of the shells. The 
relative normalizations, temperatures, and ionizations ages for each model component (mass- 
cut shell) are set from the hydrodynamics. Tying the abundances of all components together, 
we then have a model which can be fit to measured spectra by adjusting the abundances 
and an overall normalization. By using the hydrodynamics to set the plasma properties we 
reduce the degeneracy in fitting that arises when the temperature(s) and ionization age(s) 
of the plasma are unconstrained and therefore we get accurate abundance values. 
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The early deep Chandra observation, ACIS-00, is the best data set to use to determine 
the abundances of the Hll region (as opposed to the dense ER abundances, which may 
differ). These data also provide a useful example of how the fit abundances and their 
confidence ranges vary depending on the fcT-r model assumptions. The results of fitting 
four different cases are shown in Table [51 In all cases the redshifts and Nh were frozen and 
the normalization(s) and abundances (N, O, Ne, Mg, Si, S, and Fe) were fit. In the first 
case a single vpshock comp onent is fit g i ving k T ^ 3.0 keV and r ^ 1.0x10^^ scm~^; these 
are similar to the values in iPark et al.l (120021 ). As a second case, w e fit a 2-shock mode l 
fixing the temperatures, ionization ages, and Nh based on values from lZhekov et al.l ( 120101 ): 
this has a high temperature component similar to the 1-shock value along with a very low 
temperature component. In the third case, our model of Table [2] is shown, and because we 
have set A^mid to for this early data it is also a 2-shock model with a different set of fcT's 
and r's. Finally, we fit the data with the temperatures and ionization ages fixed based on 
the values from the "no ER" hydrodynamic model at the ACIS-00 epoch; these values are 
plotted in the left panel of Figure [11] and correspond to the model at the epoch shown in 
the top panel of Figure [H 

The reduced values are near 1.0 for all four of the fits, so they are equally acceptable 
from a data-fitting perspective. However, as the table shows, the l-cx ranges of the fit 
abundances depend very much on the kT-r values used in the model. If it is possible to 
constrain the kT-T values through other means or assumptions, as in the fourth case here, 
then the fit abundances will be both more realistic and better constrained. Hence, we adopt 
the hydrodynamics-based model as the most physical and include its abundances in the Hll 
column of Table [6] 
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Table 5. Degeneracy of Hll Abundances when Fitting ACIS-00 





fcT(s) 


r(s) 




Relative abundances ^ 


(1-(T range) 




Model 


(keV) 


(10" scm-3) 


N 


O 


Ne 


Si 


Fe 


1-shock 


2.94-3.34 


0.94-1.32 


0.49-0.94 


0.14-0.23 


0.00-0.17 


0.28-0.41 


0.13-0.25 


2-shock 


0.22, 2.8 


130, 1.4 


0.34-0.72 


0.11-0.17 


0.27-0.35 


0.29-0.38 


0.083-0.12 


Table |2] 


0.73, 4.2 


11.6, 5.8 


0.89-2.35 


0.28-0.63 


0.46-0.92 


0.59-0.94 


0.087-0.18 


Hydro-noER 


14 pairs 


in Figure [TT] 


0.51-0.66 


0.137-0.155 


0.35-0.42 


0.32-0.43 


0.061-0.080 



''Abundances are relative to the solar photospheric values of I Anders &: Grevessj lll989h . Abundances of Mg 
and S (not tabulated) are similar to those of Ne and Si, respectively. has been fixed at 1.3x10^""^ cm~^, unless 
stated otherwise. 

''The kT and r values for this 2-shock model were estimated from Figures 3 &; 4 of lZhekov et al. I l|20ld'l . which 
also gives an Njj of 2.23x10^"'^ cm~^ that was used in the fitting here. 



Table 6. Hydrodynamics-based Abundances and Others 





— hi: — 


— ER'* — 




Other Abundances 










Parameters 


ACIS-00 


HETG- 07&11 


Zhekov('09) 


Zhekov('lO)' 


Mattila('lO) = 


N 


0.584 ±0.07 


3.3±1.2 


-0.06 


0.32 ±0.04 


2.5±1.0 


O 


144 +0-011 
yj-i^'i'i _o.oo7 


0.382 ±0.08 


0.081 ±0.01 


0.07 ±0.01 


0.22 ±0.04 


Ne 


0.381 ±0.03 


0.566 ±0.03 


0.29 ±0.02 


0.22 ±0.01 


0.69 


Mg 


0.270 ±0.04 


0.394 ±0.02 


n no +0.01 
'J-^^ -0.02 


0.17 ±0.01 




Si 


0.373 ±0.05 


0.469 ±0.02 


n +0.02 

"■■^■^ -0.01 


0.26 ±0.01 




s 


0-293 


0.556 ±0.08 


0.30 ±0.06 


0.36 ±0.02 


0.81 ±0.4 


Ar 


0.00 +00^ 


1.2 ±0.5 


[0.537] 


[0.537] 


0.47 


Fe 


n 070 +°-°'^° 

U.U/U _o.oo9 


0.225 ±0.02 


'-'■1^ -0.00 


0.10 ±0.01 


0.20 ±0.1 


fcl u mp 




0.30 ±0.02 








Vcl u mp 




5.51 ±0.45 









Note. — Abundances are relative to the solar photospheric values of lAnders &: Grevessj l ll989ll . 
AG89. Values in brackets were fixed; those not tabulated were set to: H=l, He=2.57, C=0.03, 
Ar=0.03, Ca=0.03, and Ni=0.07. Nh has been fixed at 1.3x10^' cm-^. 

^Values are from lZhekov et all l|2009l ') Table 1 for fits to LETG and HETG data. 

''Values are from IZhekov et al ] l|2O10h Table 1 for fits to 14 ACIS spectra; a value of 
Njj=2.23±0.05 xlO^' cm~-^ was determined and used. 

'^Values are from lMattila et al. I ll20ld '). Table 7, converted to AG89 equivalents. Other of their 
values are: He=1.74 ±0.61, C undefined, and Ca=1.41. 

''The error ranges for the ER abundances and clumping parameters include the difference be- 
tween the HETG-07 and HETG-11 best-fit values and their l-cr ranges. The N abundance was set 
based on fitting the RGS-11 data over the l-cr range for ??clump- 
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Fig. 12. — Distributions of kT^-r for the hydrodynamics at the HETG-11 epoch. Points 
of {T,kTe) are plotted for all shocked mass shells in the Hll (no-ER, left panel) and "with 
ER" (right panel) simulations. Symbols for the shocked ejecta (black) and CSM (blue) are 
larger/bold for brighter shells. The Hll distribution is similar to earlier times (Figure ITTl left 
panel) although shifted to somewhat higher r and kT^ values. The "with ER" distribution 
(right panel) is dominated by components at much higher r values; the addition of a clumped 
component introduces additional /cTg-r values (red) that are scaled from the non-clumped 
CSM values. 



4.4. ER Abundances and Clumping 

We can similarly fit hydrodynamics-based models to the later observations where the 
emission is dominated by the shocked ER material and in this way determine the ER abun- 
dances. At these times there will still be a contribution from the Hll region, i.e., the v-b 
component. This contribution is removed from the data before fitting by evaluating our "no 
ER" model at the appropriate epoch and subtracting the predicted counts from the data 
before proceeding with ER-model fitting. 

The kT(.-T values of the shocked mass-cuts in our Hll and ER models at the HETG-11 
epoch are shown in Figure [121 As the blue points in the "with ER" (right) panel show, the 
kTe values in the shocked dense ring are mostly at / above 3 keV. Given these temperatures, 
simply adjusting the abundances is insufficient to fit the data, there is simply too little low- 
fcTe, high- r plasma. We have ther efore included "clumps" in the ER material following the 
scheme in iDwarkadas et al.l (l201Cl l . 



For each of the shocked CSM mass-cuts (blue circles in the right panel of Figure [T2|) we 
add an additional model component with its temperature reduced by the clump's density 
enhancement, kT^ = fcTg/r/ciump, and its ionization age increased to r' = rjciumpT- The 
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additional clumped components are shown as red circles in the right panel of Figure [12] for 
^ciump ~ 5.5. The original and clumped components have additional factors in their norms 
of (1 - /clump) and /clump ^?ciump' respectively, where /dump ~ 30% is the volume fraction of 
clumped material. The resulting ISIS fit function at the HETG-11 epoch consists of 45 NEI 
model components but with only 9 free parameters: an overall normalization, the 2 clumping 
parameters, and the ER abundances of O, Ne, Mg, Si, S, and Fe. 

We fit both the deep HETG-07 and HETG-11 data sets with their appropriate-epoch 
hydrodynamics-based models, determining the abundances and clumping parameters for 
each. There are only small differences between the 2007 and 2011 fits, so we have adopted 
nominal ER parameters that are the average values given in the "ER" column of Table [61 
Note that the N abundance was set based on similar fitting to the RGS-11 spectrum. These 
values were then used to generate our light curve (§ 14. 2p and to fit the hydrodynamics-based 
models to all of the grating spectra with only a single normalization adjustment. As exam- 
ples. Figure [13] compares the two deep HETG spectra with the sum of the hydrodynamics- 
based emission from the H II region and the nominal-parameter ER emission. 

Given our crude clumping implementation and the system's physical complexity, we 
have less confidence in the accuracy, i.e., agreement with reality, of our ER abundances 
than those we determined for the Hll region at the ACIS-00 epoch. However since the ER 
kTs and rs are based on a hydrodynamic solution there is some expectation that they will 
be an improvement over simpler mu lti-shock mod e ls. Fo r comparison, we've also included 
in Table M th e abundances given by IZhekov et al. J2OO9L multi-shock fits t o grating data), 
Zhekov et all J201ol . RSS fits to the ACIS data), and iMattila et all J201ol . optical & NIR 
data before the ER collision). Our ER abundances for N , O, Ne, and S do appear to be in 
someshat better agreement with the iMattila et al.l (|2010[ ) values than are the other X-ray- 
determined abundances. 



This is a good place to compare the free parameters of the IZhekov et al.l (120 lOl ) model 
(ZIO in the following) with ours, in terms of their number and their values. The two model 
geometries are similar at later times, both are dominated by the ER emission. Figure [6] 
and both make approximations to include clumped ER material. Each model fixes time- 
independent values for the CSM's p(r), the N^, and the abundances. In addition to these, 
ZIO fit six parameters for each observation: the values of the temperature, ionization age, 
and emission measure for both the blast wave (BW) and the transmitted shock (TS). For 
example, for their latest data point which is near the HETG-07 epoch, the values for these 
are: kT^w -1-7 keV, kT^s -0.33 keV, tbw -1.4x10^1 scm^^^ tbw -30x10" scm'^, and 
the respective emission measures are 1.05 and 9.5 in units of lO^^cm"'^. In contrast, all of 
our mass-shell kT and r values are fixed by the hydrodynamics, e.g., as for HETG-11 in 
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Figure [T2|, and our ER modeling adds only three further time-independent parameters: the 
overall normalization, Qer, and the two clumping parameters. In our approximation the 
clump over-density, r^ciump ~ 5.5, sets the ratio of ER-to-clump temperature. For the ZIO 
values above, this ratio is Tbw/^ts ~ 5, close to ours. Likewise, the ZIO emission measure 
ratio, EMts/EMbw ~ 9, is in the ball park of our value which is given by /ciump^ciump/(l ~ 
/clump) ~ 13, using our clump volume fraction /dump ~ 30 %. In constrast, the clump volume 
fraction for the ZIO model is of order 5 % because the clump-to- ambient density ratio is large: 
Pts/pbw ~ 14; t his latter valu e is a f unction of the temperature ratio as shown in a plot in 
the Appendix of IZhekov et al.l (120091 ). These comparisons suggest that there may be some 
degeneracy between over-density and volume fraction for the clumps. In anycase, it is clear 
that both approaches here are only approximating the complex 2D/3D density structures of 
the real ER. 
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Fig. 13. — HETG(MEG) spectra and full hydrodynamics-based models. The HETG-07 
(top) and HETG-11 (bottom) data (black) are reasonably fit by the total (Hll plus nominal 
clumped-ER) model spectra (red). The HlI-shocked-CSM (blue) and the Hll-shocked-ejecta 
(green) components are also shown individually. Note the relative growth of the shocked 
ejecta (green) between the epochs; hence at late times the ejecta and its abundances will 
dominate the Hll (i.e., v-b) component. 



- 39 - 




1 2 
Energy [keV] 



Fig. 14. — Comparing the heuristic very-broad spectrum with the hydrodynamics-based Hll 
v-b spectrum. The HETG-07 data (black) are shown with the best-fit 3-shock model (red) 
that includes a v-b component (blue). This "smoothed version of the model" form of the 
v-b component is compared with the multi-shells model spectrum at the HETG-07 epoch 
(magenta). Ignoring the overall continuum shape, they are in reasonable line-broadening 
agreement except for the lack of significant FeXVII lines between 0.70-0.75 keV in the hydro- 
based model. 

4.5. Very-Broad Spectrum from the Simulation 

We can check our decision to use a smoothed version of the spectrum as the model for 
the v-b component by directly comparing the two of them, as shown in Figure [TH Here, the 
broad lines in the "real" (from the hydrodynamics) v-b component are similar to those in 
the smoothed version of the 3-shock fit to the data. This similarity can be roughly explained 
using the kT^-r distributions shown in Figure the H II (v-b) emission is dominated by 
components with kT^ ~ 2.5 keV and r ^ GxlO^*^ scm~^, whereas the main ER emission is 
from clumped material with kT^ ~ 0.7 keV and r ^ 2x10^^ scm~^. Putting these values 
into a vnei model gives similar line structures for the two cases, hence, the high-/cT-low-r 
plasma produces ionization states and emission lines that are similar to those of a low-fcT- 
high-r plasma. Based on this comparison we did, however, decide to ignore the FeXVII 
range 0.70-0.75 keV when doing the v-b fitting as these lines were clearly not present in the 
hydro dynamic-based v-b spectrum. 

In terms of the assumption of a constant width of the v-b emission in time (the evolution 
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path shown in Figure H]), the average bulk motion of the shocked Hll material in the "no 
ER" simulation varies (only) from 4200 kms~^ to 3850 kms~^ from the HETG-99 epoch to 
an age of over 25 years (9000+ days). This is also seen in the locations of the "no ER" CD 
and FS locations, lower-left panel of Figure 13 which follow nearly straight lines with overall 
slopes of 4300 kms~^ and 5100 kms~^, respectively. 



4.6. The Mass of the X-Ray Emitting Material 

We can calculate the shocked mass for different components and epochs of the model. 
For the Hll region, at the ACIS-00 epoch the FS is at ~ 0.162 pc (top panel of Figure [H]) and 
the mass of shocked CSM in the simulation, correcting for the opening angle normalization 
OHii/47r ^ 0.26 (§ 14. 2p . is 0.012 solar masses. At this time our model shows ~ 0.005 solar 
masses of ejecta are shocked within the same opening angle. Later, at the HETG-11 epoch, 
the shocked Hll material that produces the v-b emission in our model has grown to 0.042 
solar masses along with an additional 0.024 solar masses of shocked ejecta. Our Hll mass 



values compare in order of magnitude with the mass of 0.018 solar masses that iMattila et al. 



(120 lOl ) require in their low-density (10^ atoms cm ^) Hll component. 



The total "47r" mass of the thin ER is 1.18 solar masses; however, correcting for the 
small solid-angle fraction, Q-eb./4:Tt ~ 0.0016, we get only ~ 0.0019 solar masses of ER that 
are shocked at the HETG-11 epoch (when the FS is exiting the thin ER). Including the extra 
mass due to clumping (multiplying by a factor of 2.35, § 14.21) . we have 0.0045 solar masses 
as the total shocked-ER mass that produces the observed non-v-b X-rays in our model at 
the HETG-11 epoch. This is to be compared with an estimate of the total mass of the 
UV-ionized ER of 0.058 solar masses distrib uted among three densities of 1.4x10^, 4.2x10^, 



and 4.2x10"^ amucm"^ (IMattila et al.ll2010l . using ~ 1.4 amu per atom); of this total, 0.0120 
solar masses are in the highest-density component, almost 3 times our X-ray inferred mass. 
However, it could well be that our ER X-ray emission is a factor of a few overly efficient per 
mass and including this as an error bar puts us in a gray area in trying to use this mass 
comparison to decide if all of the ionization-visible dense ER has been shocked at the current 
~ 25 year age of SN 1987A. 



4.7. Image Implications 



We can construct schematic geometric visualizations of SN 1987A using our simple hy- 
drodynamic models to set the locations and fluxes of the components at a given epoch. The 
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main emission components at two different epochs are shown in multiple views in Figure [T5] 
using a color-coding similar to that for the shocked plasma in th e schematic of Figure l6] 
These images were created using simple software extensions to ISIS ( Dewey fc Noblell2009l )rl. 



The Hll ("no ER") simulation gives rise to emission from i) shocked Hll material and 
ii) reverse- shocked ejecta; these regions extend out of the ER plane as shown in the images. 
The "with ER" simulation is shown using a schematic geometry of 21 discrete spots to 
emphasize our "X-ray hot spots" conclusion that is based on the small value of ^^er- The 
"with ER" emission components include: iii) the shocked dense ER and its clumps, iv) the 
shocked Hll material that is in the plane of the ER, and v) reverse- shocked ejecta also in the 
ER plane. Because the effective solid angle associated with the ER simulation is so small 
(^ER/4vr = 0.0016) these last two components make only minor flux contributions which are 
exceeded by the first three. We note that component "iv)" is the plasma that is "shocked 
again" in the reflected shock structure (RSS) paradigm (IZhekov et al.ll2009U2010l ): our "with 
ER" hydrodynamics does show the reflected shock(s) having velocity effects on the shocked 
HlI region interior to the shocked ER for about 1.5 years, after which the region's velocity 
recovers and is similar to other regions, e.g., see the velocity curve in Figure |S1 bottom panel. 

These visualizations suggest the need for fitting the X-ray imaging data with multiple 
sp atial coniponent s and provide some guidance on the component properties. Recent work 
of iNg et al.l (|2009[ ) has explored a two-component spatial model in fitting SN 1987A's image, 
and as they emphasize, further observations are needed to come to a firm conclusion. 



See also: |http : //space . mit . edu/home/dd/Event2D/ 
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Fig. 15. — Hydrodynamics-based geometric emission models of SN 1987A at epochs 2004.5 
(left column) and 2011.5 (right column). Three views are shown: a "polar" view (top; 
looking perpendicular to the equatorial plane), a "side" view (middle; north is up and the 
Earth is to the left), and the usual sky view (bottom; north up, east to left). Four emission 
components are color-coded: the shocked Hll material (blue), the reverse-shocked ejecta 
(green), the emission from the shocked ER (red spots), and higher-lattitude 8.6 GHz radio 
emission (gray). 
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5. CONCLUSIONS 

In this paper we have modelled the X-ray emission from SN 1987A, focussing on the 
very-broad (v-b) component seen with the high-resolution grating instruments. Although 
the CSM of SN 1987A is known to have a complex shape from optical observations, we have 
chosen to model it using ID simulations in and out of the equatorial plane, our Hll and 
ER hydrodynamics. Even with these simplifying assumptions the models yield a good fit 
to the X-ray-measured radii and growth rates, the multi-band X-ray light curves, and the 
high-resolution X-ray spectra. 

Our results show that better fits to the X-ray spectra over the last decade are obtained 
if a very-broad component of emission, with a FWHM of about 9300 km s'S is included. 
Such a v-b component is present in the data over at least the last decade. Our hydrodynamic 
simulations provide a natural explanation for this component: it arises from the shocked Hll 
material and it is the dominant X-ray component until ~ 5500 days after the SN explosion. 
Since then the total 0.5-2 keV flux dramatically increased due to the "ER" collision, yet the 
inferred 0.5-2 keV flux of the v-b component itself continued to grow in agreement with our 
hydrodynamics expectations. Identifying the 3-10 keV flux as originating primarilly from the 
v-b component provides a natural explanation for the difference in the 0.5-2 and 3-10 keV 
light curves. At the present epoch, the v-b contribution is ~20% of the total 0.5-2 keV flux, 
and of this the hydrodynamics suggests that roughly half is coming from reverse-shocked 
ejecta. This ejecta contribution should grow and may cause the v-b component to become 
selectively enhanced in lines that reflect the (outer) ejecta composition. 

Our ER hydrodynamics adquately reproduces the observed X-ray image radii and the 
bulk of the late-time 0.5-2 keV flux increase. We do require clumps in the ER with den- 
sity enhancements of 5.5 x the ~ 10^ amucm"^ ambient ER values, not unlike the clumps 



that IZhekov et al.l (120101 ) have added within their smooth CSM density profile. These are 
approximations to the actual 2D/3D density structure of the shock-ER interaction; it would 
not be surprising if the effective clumping parameters might change in time, e.g., as clumps 
evaporate and/or turbulent structures evolve. 

Our approach allows us to make predictions of what the image size. X-ray spectrum and 
light curve will look like in the next few years, in particular for the two cases of i) a continuing 
on-going interaction with dense ER material (our "thick ER" case) and ii) the drop in flux 
if the ER is "thin" . In this latter case, the FS in the equatorial plane has generally exited 
th e densest pa r ts of the ER and has begun moving into lower-density material as suggested 



by iPark et al.l ( 1201 ll ). For this case our hydrodynamics shows a decrease in the 0.5-2 keV 
flux of ~17% per year. Note that a third option that cannot be completely ruled out is 
that the FS will encounter larger amounts of dense ER material, e.g., the main body of the 
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ER; this would result in even greater increases in flux than we have shown here. Comparing 
these and more realistic multi-dimensional predictions with future observations will provide 
significant information regarding the thickness, density, and structure of the equatorial ring 
and the Hll region within and outside the equatorial plane as SN1987A moves into a new 
phase. 
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FaciUUes: CXO (HETG). 

A. Spectral Models and the gsmooth Function 

The 3-shock model used in § l2.3l makes use of the XSPEC spectral library and is defined 
in ISIS with the command: 

f it_fun( "phabs(l) * gsmooth(l, vpshock(l) + vpshock(2) + vpshock(3) )" ); (Al) 

where the trailing ";" is the S-Lang end-of-expression character. The gsmooth (1) function 
here is used to add blur to the model in order to account for spatial and spectral effects that 
are unique to SN 1987A; the gsmooth parameters Sig@6keV and Index and their values for 
this purpose are described further below. 

A similar expression is used to define the v-b 3-shock model used in § 12.41 as the sum of 
the original 3-shock model and a very-broad version of itself: 



fit_fun( "setvb(l) * phabs(l) * " -f- (A2) 
" ( (1 — f racbroad) * gsmooth(l, vpshock(l) + vpshock(2) + vpshock(3)) -|- " -|- 
"f racbroad * gsmooth(2, vpshock(l) + vpshock(2) + vpshock(3)) )" ); 
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Here the S-Lang variable, fracbroad, sets the fraction of the total which is in the v-b 
component and the second gsmooth term sets the width of the v-b emission. For a v-b 
component having a FWHM of f^.b, the Sig@6keV parameter is set to 6fv-b/(2.35c) and the 
Index parameter is set to 1, as discussed further below. The model component (setvb) is a 
dummy function which returns unity in the spectral model but is used to control the broad- 
fraction variable from an ISIS fit parameter, setvb (1) .fb. This is accomplished before the 
model is defined with a few lines of code: 

public variable fracbroad; 

% Define a dummy model function to set the broad fraction 
define setvb_fit (l,h,p) 

{ 

% set fracbroad from the input parameter: 
fracbroad = p [0] ; 
return 1; 

} 

add_slang_function ("setvb", ["fb"]); 

In the remainder of this section we provide more specifics on the gsmooth function and its 
parameters. 

The gsmooth function convolves a spectral model with a Gaussian function whose width, 
a El varies with energy according to: 

(Je = (Te (^kev/6.0)" . (A3) 

The constants of this equation, ue and a, are set through the values of the gsmooth param- 
eters Sig@6keV and Index, respectively. This versatile function can be used to approximate 
instrumental effects in place of an rmf . For the case of a = it produces a const ant- with- 
energy blur, e.g., simulating a micro-calorimeter-like response. The case of a = 0.5 gives a 
resolution, (i£^FWHM, that goes as the square-root-of-energy, e.g., as in a proportional-counter 



Table 7. Spatial-Spectral Blur Parameters 





CO VIII 


o-Mg 




aa, SigaekeV 


Grating 


(eV) 


(cV) 


Q, Index 


(keV) 


RGS 


0.628 


[1.36] 


1.0000 


0.00580 


HETG 


0.881 


2.14 


1.1454 


0.01124 


LETG 


1.893 


5.90 


1.4682 


0.04946 
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or a CCD detector. For a = 1 the resolving-power, E/dE-pM^mAi is constant, and for a = 2 
we have the case of a blur that is constant in wavelength as is nearly the case for a grating 
instrument like the HETG. 

The smoothing term can also be used to include additional source spectral properties, for 
example, the Doppler blur from a Gaussian velocity spread, frms, is obtained by setting a = 1 
and o"6 = 6 frms/ c ; if the FWHM velocity width is given, then we have Ue = 6 frwHivi/ (2.35c). 
It is in exactly this way that we use gsmooth(2) in Expression IA2I to produce the v-b 
component in our model. 

Finally, we also make use of the gsmooth function to account for additional blurs that 
are present in the narrow lines of the SN 1987A grating spectra but are not included in 
the spectral model (vpshock) nor in the point-source instrumental resp onses (rmfs). These 



blurs have been rioted and included in previous grating-data analyses (IZhekov et al.l 12005 



Dewey et al.ll2008l : ISturm et al.ll2010l ) and result from i) bulk and turbulent motion of the X- 
ray emitting plasma and ii) spatial-spectral effects which broaden the observed line widths 
due to the spatial size of SN 1987A. In the context of the gsmooth function the former 
is accounted for with an a = 1 behavior and the latter with an a = 2 dependence. The 
combination of both effects can be approximated using an intermediate value of a determined 
by measuring the narrow-component line width at two energies. For the analyses here this 
has been done using regions including the OVIII line and the Mg line (1.28 - 1.55 keV); the 
measured values and resulting gsmooth parameters are given in Table [71 Note that the RGS 
has a = 1 because SN 1987A is spatially unresolved by XMM-Newton. 
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